!
! astrofits - write WCS to FITS image & read star table
!
!
! Copyright © 2011 F.Hroch (hroch@physics.muni.cz)
!
! This file is part of Munipack.
!
! Munipack is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
! 
! Munipack is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
! 
! You should have received a copy of the GNU General Public License
! along with Munipack.  If not, see <http://www.gnu.org/licenses/>.
!

module astrofits

  use fitsio

  implicit none

  integer, parameter, private :: dbl = selected_real_kind(15)

contains

  subroutine readtable(filename, key_alpha, key_delta, alpha, delta, rframe, catid)

    character(len=*), intent(in) :: filename, key_alpha, key_delta
    real(dbl), allocatable, dimension(:), intent(out) :: alpha, delta
    character(len=*), intent(out) :: rframe, catid

    integer :: status, nrows
    logical :: anyf
    integer :: frow = 1, felem = 1
    real(dbl) :: nullval = 0.0_dbl
    character(len=80) :: buf, colname1, colname2
    integer :: col1, col2

    status = 0

    ! open and move to a table extension
    call fttopn(15,filename,0,status)
    call ftgnrw(15,nrows,status)
    if( status /= 0 ) goto 666

    ! define reference frame and identification of catalogue
    rframe = ''
    catid = ''

    call ftgcnn(15,.false.,key_alpha,colname1,col1,status)
    call ftgcnn(15,.false.,key_delta,colname2,col2,status)
    if( status /= 0 ) goto 666

    if( nrows > 0 ) then
       allocate(alpha(nrows),delta(nrows))
       call ftgcvd(15,col1,frow,felem,nrows,nullval,alpha,anyf,status)
       call ftgcvd(15,col2,frow,felem,nrows,nullval,delta,anyf,status)
    end if
    call ftclos(15,status)

    if( .not. (nrows > 0) ) then
       write(*,*) 'Empty catalogue table.'
       stop 1
    end if

666 continue
    if( status /= 0 )then
       write(*,*) 'FITSIO: following error(s) has been occured:'
       do
          Call ftgmsg(buf)
          if( buf == ' ' ) exit
          write(*,*) trim(buf)
       enddo
       stop 1
    endif

  end subroutine readtable


  subroutine readfile(filename, w, h, x, y, crpix)

    character(len=*), intent(in) :: filename
    real(dbl), intent(out) :: w, h
    real(dbl), allocatable, dimension(:), intent(out) :: x, y
    real(dbl), dimension(:), optional, intent(out) :: crpix

    integer :: status, blocksize, naxis, pcount,gcount, bitpix, nrows
    integer, dimension(2) :: naxes
    logical :: simple,extend,anyf
    integer :: extver = 0, frow = 1, felem = 1
    real(dbl) :: nullval = 0.0_dbl
    character(len=80) :: buf, colname1, colname2
    integer :: col1, col2

    status = 0
    call ftopen(15,filename,0,blocksize,status)
    if( status /= 0 ) goto 666 !stop 'Failed to open FITS file.'

    call ftghpr(15,2,simple,bitpix,naxis,naxes,pcount,gcount,extend,status)
    if( status /= 0 .or. naxis /= 2 ) goto 666
    w = naxes(1)
    h = naxes(2)
    if( present(crpix) ) then
       crpix = naxes/2.0_dbl
    end if

    call ftmnhd(15,BINARY_TBL,MEXTNAME,extver,status)
    if( status == BAD_HDU_NUM ) then
       write(*,*) 'MUNIPACK extension not found. A photometry is not available yet.'
       stop 1
    end if
    call ftgnrw(15,nrows,status)
    if( status /= 0 ) goto 666
    allocate(x(nrows),y(nrows))

    ! the column labels are specific for MUNIPACK extension only
    call ftgcnn(15,.false.,'x',colname1,col1,status)
    call ftgcnn(15,.false.,'y',colname2,col2,status)

    call ftgcvd(15,col1,frow,felem,size(x),nullval,x,anyf,status)
    call ftgcvd(15,col2,frow,felem,size(y),nullval,y,anyf,status)
    call ftclos(15,status)

666 continue
    if( status /= 0 )then
       write(*,*) 'FITSIO: following error(s) has been occured:'
       do
          Call ftgmsg(buf)
          if( buf == ' ' ) exit
          write(*,*) trim(buf)
       enddo
       stop 1
    endif
    
  end subroutine readfile



!-----------------------------------------------------------------------------

  subroutine wcsupdate(file,output,type,crval,dcrval,sc,dsc,pa,dpa,fcrpix,crpixels,com)

    integer, parameter :: ndim = 2
    real(dbl), parameter :: rad = 57.295779513082322865_dbl

    character(len=*), intent(in) :: file, output, type
    real(dbl), dimension(:), intent(in) :: crpixels,crval,dcrval
    real(dbl), intent(in) :: sc,dsc,pa,dpa
    logical, intent(in) :: fcrpix
    character(len=*), dimension(:), intent(in), optional :: com

    integer :: status, blocksize, naxis, pcount,gcount, bitpix
    integer, dimension(ndim) :: naxes
    real(dbl), dimension(ndim) :: crpix
    logical :: simple,extend
    real(dbl) :: a,b,da,db
    character(len=80) :: buf
    integer :: i

    a = cos(pa/rad)
    da = (dpa/rad)*sin(pa/rad)
    b = sin(pa/rad)
    db = (dpa/rad)*cos(pa/rad)

    status = 0
    
    if( file /= '' .and. output /= '' .and. file /= output ) then

       ! backup
       call ftopen(15,file,1,blocksize,status)
       call ftinit(16,output,blocksize,status)
       call ftcpfl(15,16,1,1,1,status)
       call ftclos(16,status)
     
    else if( output /= '' .and. output == file ) then
     
       ! overwrite
       call ftopen(15,output,1,blocksize,status)

    else

       write(*,*) 'Internal file handle error.'
       stop 666

    end if

    if( status /= 0 ) goto 666

    if( fcrpix ) then
       crpix = crpixels

    else

       call ftghpr(15,2,simple,bitpix,naxis,naxes,pcount,gcount,extend,status)
       if( status /= 0 ) goto 666
       crpix = naxes/2.0_dbl
       
    end if

    if( type == "GNOMONIC" )then
       call FTUKYS(15,'CTYPE1','RA---TAN','the coordinate type for the first axis',status)
       call FTUKYS(15,'CTYPE2','DEC--TAN','the coordinate type for the second axis',status)
    end if
    call FTUKYD(15,'CRVAL1',crval(1),16, &
         '[deg] first axis value at reference pixel',status)
    call FTUKYD(15,'CRVAL2',crval(2),16, &
         '[deg] second axis value at reference pixel',status)
    call FTUKYD(15,'CRPIX1',crpix(1),16, &
         'x-coordinate of reference pixel',status)
    call FTUKYD(15,'CRPIX2',crpix(2),16, &
         'y-coordinate of reference pixel',status)
    call FTUKYD(15,'CD1_1',-a/sc,16, &
         'partial of first axis coordinate w.r.t. x',status)
    call FTUKYD(15,'CD1_2',b/sc,16, &
         'partial of first axis coordinate w.r.t. y',status)
    call FTUKYD(15,'CD2_1',b/sc,16, &
         'partial of second axis coordinate w.r.t. x',status)
    call FTUKYD(15,'CD2_2',a/sc,16, &
         'partial of second axis coordinate w.r.t. y',status)

    if( present(com) ) then
       do i = 1, size(com)
          call ftpcom(15,com(i),status)
       end do
    end if

    call ftclos(15,status)

666 continue

    if( status /= 0 )then
       write(*,*) 'FITSIO: following error(s) has been occured:'
       do
          Call ftgmsg(buf)
          if( buf == ' ' ) exit
          write(*,*) trim(buf)
       enddo
       stop 1
    endif

  end subroutine wcsupdate

end module astrofits
